library(xtable)

ploughs <- foreign::read.dta("crosscountry_dataset.dta")

ploughs$centered_ln_inc <- ploughs$ln_income - mean(ploughs$ln_income, na.rm = TRUE)
ploughs$centered_ln_incsq <- ploughs$centered_ln_inc^2

ate.mod <- lm(women_politics ~ plow + agricultural_suitability +  tropical_climate +  large_animals + political_hierarchies + economic_complexity + rugged, data = ploughs)
summary(ate.mod)

naive.mod <- lm(women_politics ~ plow * centered_ln_inc + plow * centered_ln_incsq + agricultural_suitability +  tropical_climate +  large_animals + political_hierarchies + economic_complexity + rugged + centered_ln_inc + centered_ln_incsq, data = ploughs)
summary(naive.mod)

ptbias.mod <- lm(women_politics~ plow*centered_ln_inc + plow*centered_ln_incsq + agricultural_suitability +  tropical_climate +  large_animals + political_hierarchies + economic_complexity + rugged +years_civil_conflict +  years_interstate_conflict  + oil_pc + european_descent + communist_dummy +polity2_2000+serv_va_gdp2000, data = ploughs)
summary(ptbias.mod)

med.mod <- lm(centered_ln_inc ~ plow + agricultural_suitability +  tropical_climate +  large_animals + political_hierarchies + economic_complexity + rugged, data = ploughs)
summary(med.mod)


direct.mod <- lm(I(women_politics - coef(ptbias.mod)["centered_ln_inc"]*centered_ln_inc - coef(ptbias.mod)["centered_ln_incsq"]*centered_ln_incsq - coef(ptbias.mod)["plow:centered_ln_inc"]*plow*centered_ln_inc - coef(ptbias.mod)["plow:centered_ln_incsq"]*plow*centered_ln_incsq) ~ plow + agricultural_suitability +  tropical_climate +  large_animals + political_hierarchies + economic_complexity + rugged, data = ploughs)
summary(direct.mod)

nboots <- 1000
ate.hold <- rep(NA, times = nboots)
naive.hold <- rep(NA, times = nboots)
ptbias.hold <- rep(NA, times = nboots)
acde.hold <- rep(NA, times = nboots)
set.seed(02143)
for (i in 1:nboots) {
    star <- sample(1:nrow(ploughs), replace = TRUE)
    ploughs.star <- ploughs[star,]

    ate.star <- lm(women_politics ~ plow + agricultural_suitability +  tropical_climate +  large_animals + political_hierarchies + economic_complexity + rugged, data = ploughs.star)
    naive.star <- lm(women_politics ~ plow*centered_ln_inc + plow*centered_ln_incsq + agricultural_suitability +  tropical_climate +  large_animals + political_hierarchies + economic_complexity + rugged + centered_ln_inc + centered_ln_incsq, data = ploughs.star)
    ptbias.star <- lm(women_politics~ plow*centered_ln_inc + plow*centered_ln_incsq + agricultural_suitability +  tropical_climate +  large_animals + political_hierarchies + economic_complexity + rugged +years_civil_conflict +  years_interstate_conflict + oil_pc + european_descent + communist_dummy +polity2_2000+serv_va_gdp2000, data = ploughs.star)
    direct.star <- lm(I(women_politics - coef(ptbias.star)["centered_ln_inc"]*centered_ln_inc - coef(ptbias.star)["centered_ln_incsq"]*centered_ln_incsq - coef(ptbias.star)["plow:centered_ln_inc"]*plow*centered_ln_inc - coef(ptbias.star)["plow:centered_ln_incsq"]*plow*centered_ln_incsq) ~ plow + agricultural_suitability +  tropical_climate +  large_animals + political_hierarchies + economic_complexity + rugged, data = ploughs.star)

    ate.hold[i] <- coef(ate.star)["plow"]
    naive.hold[i] <- coef(naive.star)["plow"]
    ptbias.hold[i] <- coef(ptbias.star)["plow"]
    acde.hold[i] <- coef(direct.star)["plow"]
}

cat("ATE: ", mean(ate.hold), " [", quantile(ate.hold, 0.025), ", ", quantile(ate.hold, 0.975), "]\n",
    "Naive: ", mean(naive.hold), " [", quantile(naive.hold, 0.025), ", ", quantile(naive.hold, 0.975), "]\n",
    "ACDE: ", mean(acde.hold), " [", quantile(acde.hold, 0.025), ", ", quantile(acde.hold, 0.975), "]\n")


acde <- coef(direct.mod)["plow"]
acde.ci <- quantile(acde.hold, c(0.025, 0.975))

ate <- coef(ate.mod)["plow"]
ate.ci <- quantile(ate.hold, c(0.025, 0.975))

naive <- coef(naive.mod)["plow"]
naive.ci <- quantile(naive.hold, c(0.025, 0.975))


cis <- sprintf("%.2f", c(ate.ci, naive.ci, acde.ci))
cis <- apply(matrix(cis, ncol = 2, byrow = TRUE), 1, paste, collapse = ", ")
cis <- paste("[", cis, "]", sep = "")
res <- data.frame(rbind(ate, naive, acde), cis)
rownames(res) <- c("ATE", "Ignoring $Z_i$", "Sequential g")
colnames(res) <- c("Estimate", "95\\% Bootstrapped CI")
xres <- xtable(res, align = "lll")
print(xres, floating = FALSE, sanitize.text.function=identity)


biased.diff <- ate.hold - naive.hold
adie <- ate.hold - acde.hold


#cairo_pdf(file = "../../figs/ploughs.pdf",width = 4.5, height = 3.5, pointsize = 9, family = "Minion Pro")
plot(density(adie), xlim = c(-7,14), ylim = range(density(biased.diff)$y), lwd = 2, yaxt = "n", main = "Strength of Mechanism (ATE-ACDE)", xlab = "Percent Political Positions Held by Women", bty = "n", ylab = "Bootstrap distribution")
lines(density(biased.diff), col = "grey70", lwd = 2)
axis(side = 2, las = 2)
abline(v = 0, lty = 2, col = "grey30")
text(x = -3.75, y = 0.076, "Ignoring Intermediate\nConfounders", col = "grey70")
text(x = 12.5, y = 0.076, "Sequential\ng-estimation")
#dev.off()


nboots <- 1000
incvals <- seq(4.5, 10.5, by = 0.5)
acde.bmat <- matrix(NA, nrow = nboots, ncol = length(incvals))
set.seed(02143)
for (i in 1:nboots) {
  star <- sample(1:nrow(ploughs), replace = TRUE)
  ploughs.star <- ploughs[star,]
  ate.star <- lm(women_politics ~ plow + agricultural_suitability +  tropical_climate +  large_animals + political_hierarchies + economic_complexity + rugged, data = ploughs.star)

  for (j in 1:length(incvals)) {
    ## recenter the M variable
    ploughs.star$thisinc <- ploughs.star$ln_income - incvals[j]
    ploughs.star$thisincsq <- ploughs.star$thisinc^2
    ptbias.star <- lm(women_politics~ plow * thisinc + plow * thisincsq + agricultural_suitability +  tropical_climate +  large_animals + political_hierarchies + economic_complexity + rugged + years_civil_conflict +  years_interstate_conflict  + oil_pc + european_descent + communist_dummy +polity2_2000+serv_va_gdp2000, data = ploughs.star)
    direct.star <- lm(I(women_politics - coef(ptbias.star)["thisinc"]*thisinc - coef(ptbias.star)["thisincsq"]*thisincsq - coef(ptbias.star)["plow:thisinc"]*plow*thisinc - coef(ptbias.star)["plow:thisincsq"]*plow*thisincsq) ~ plow + agricultural_suitability +  tropical_climate +  large_animals + political_hierarchies + economic_complexity + rugged, data = ploughs.star)
    acde.bmat[i,j] <- coef(direct.star)["plow"]
  }
}


#cairo_pdf(file = "../../figs/ploughs-interaction.pdf",width = 4.5, height = 3.5, pointsize = 9, family = "Minion Pro")
plot(x = incvals, y = colMeans(acde.bmat), pch = 19, ylim = range(acde.bmat), bty = "n", yaxt = "n",
     xlab = "Log GDP per capita, 2000", ylab = "ACDE of the Plough")
axis(side = 2, las = 2)
abline(h = 0, col = "grey")
segments(x0 = incvals, y0 = apply(acde.bmat, 2, quantile, 0.025), y1 = apply(acde.bmat, 2, quantile, 0.975))
rug(x = ploughs$ln_income)
#dev.off()
